Projeto da disciplina de Eletromagnetismo

Alunos:

Kennedy Edmilson Cunha Melo kecm@cin.ufpe.br

Rafael dos Reis de Labio rrl3@cin.ufpe.br

Aplicação do Método dos Momentos para determinar a distribuição superficial de carga nas placas abaixo

# Importando as bibliotecas
import numpy as np
import matplotlib.pyplot as plt
import math

Pontos centrais de cada elemento (condição de contorno pela divisão do lado maior em N segmentos, em que N é par):

Pontos dados por: $ (x, y, z) = (i \cdot \Delta - \frac{\Delta}{2}, j \cdot \Delta - \frac{\Delta}{2}, k \cdot d) $ onde $ \Delta = \frac{L}{N} $ , $ k = 0, 1 \,$ ; $\, 1 \leq i \leq N\,$ e $\, 1 \leq j \leq N \,$ para $ i \leq \frac{N}{2} \,\, e \,\, j\leq \frac{N}{2}.\,\,$ já para $i \geq \frac{N}{2} \,\, e \,\, j\geq \frac{N}{2}$ Os pontos centrais são considerados nulos

# Variáveis iniciais

e0 = 8.85e-12   # Permissividade no vácuo
Vi = 0          # Potencial da placa em z = 0
V0 = 10         # Potencial da placa em z = d
d = 0.002       # A distancia entre as placas condutoras = 2 mm
L = 0.1         # O Lado = 10 cm
N = 20          # Número de segmentos
delta = L/N     # Tamanho do lado de cada elemento 
                # Cada lado é divido em N segmentos, totalizando N² elementos
                # Porém, não iremos consiserar 1/4 faltante da placa, com isso teremos 3*((N²)/4) de cada placa
qtdElementos = int(3*((N**2)/2)) #quantidade de discretizações das duas placas serão 3*((N²)/4) + 3*((N²)/4) = 3*((N²)/2)

 
#Estes valores das constantes serão para a questão 2, porém iremos variar o N.

# Matriz para as coordenadas do centro de cada elemento (x, y, z)
# A matriz vai ter 2 * N² elementos (duas placas condutoras)

def coordenadasCentro(n):
    matrizPosicao = []
    
    for k in range(2):
        for i in range(n):
            for j in range(n):
                # checamos se no eixo y <= L/2 para discretizarmos todas as células delta, 
                # caso não, só discretizamos com y > L/2, apenas com x <= L/2
                val = [(j * delta) + (delta/2), (i * delta) + (delta/2), k * d]
                if i >= (n/2) and j >= (n/2):
                  continue
                else:
                  matrizPosicao.append(val)
                  
    return matrizPosicao

matrizPosicao20 = coordenadasCentro(N)
#Printando matriz
# for vet in matrizPosicao20:
#   for val in vet:
#     print(round(val, 3), end = ' ')
#   print()

Elementos da matriz de impedância [Z]:

Os elementos dessa matriz são da forma $ Z_{mn} $ (onde: $ 1 \leq m \leq \frac{3N^{2}}{2} $ e $ 1 \leq n \leq \frac{3N^{2}}{2} $), e aqui temos dois casos:

Se $ m = n $:

$ Z_{mn} = \frac{\Delta}{\pi \epsilon_{0}}ln(1 + \sqrt{2}) $

Se $ m \neq n $, podemos aproximar o valor da impedância como se toda a carga estivesse concentrada no centro:

$ Z_{mn} \approx \frac{1}{4 \pi \epsilon_{0}} \frac{\Delta^{2}}{\sqrt{(x_{m} - x_{n})^{2} + (y_{m} - y_{n})^{2} + (z_{m} - z_{n})^{2}}} $

# Matriz de impedância
def impendancia(m, n, R):
    # R = matrizPosicao
    if m != n:
        # 1/(4 * pi * e0) * delta²/sqrt[(xm - xn)² + (ym - yn)² + (zm - zn)²]
        return 1/(4 * math.pi * e0) * (delta**2)/math.sqrt(((R[m][0] - R[n][0])**2) + ((R[m][1] - R[n][1])**2) + ((R[m][2] - R[n][2])**2))
    
    # [delta/(pi * e0)] * ln(1 + sqrt(2))
    return delta/(math.pi * e0) * math.log(1 + math.sqrt(2))
    
def matrizImpedancia(n, matrizPosicao):
    Z = []
    # tamanho = dimensão da placa L de baixo (Z=0) que é (3/4) * n² + dimensão da placa L de cima (z=d) que é (3/4) * n²
    for i in range(n):
        Z.append([])
        for j in range(n):
            Z[i].append(impendancia(i, j, matrizPosicao))
      
    return Z

# Matriz de impedância para o caso inicializado (N == 20)          
Z20 = matrizImpedancia(qtdElementos, matrizPosicao20)
  
#Printando a matriz
# for i in range(qtdElementos):
#     for j in range(qtdElementos):
#         print(round(Z20[i][j], 3), end = ' ')
#     print()

Elementos da matriz de tensão [V]:

$ V_{m} = 0, 0 \leq m \leq \frac{3N^{2}}{4}\, \, $ ou $\, \, V_{m} = V_{0}\,, \,\frac{3N^{2}}{4} < m \leq \frac{3N^{2}}{2}\,\,$e nos pontos os quais a placa existe, obedecendo seu formato em L

# Matriz de tensão (potencial)

def matrizTensao(n):
    V = []
    # // é divisão por inteiro
    for i in range(n//2): # Placa condutora em z = 0
        V.append(Vi)
        
    for i in range(n//2, n): # Placa condutora em z = d
        V.append(V0)
        
    return V

# Matriz de tensão para o caso inicializado (N == 20)
V20 = matrizTensao(qtdElementos)
 
# Printando a matriz
for i in range(qtdElementos):
    print(V20[i], end = ' ')

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 

Segunda questão

def amplitude(imp, volt):
    # Resolvendo o sistema para determinar as amplitudes dos pulsos
    return np.linalg.solve(imp, volt)
    # Cada elemento de A20 terá uma solução para cada parte dividida da placa condutora (que está em R)

# Amplitude dos pulsos para o caso inicializado (N == 20)
A20 = amplitude(Z20, V20)
# print(A20)

def distSuperf(amp):
    # A distribuição superficial de carga pode ser aproximada por uma soma dos elementos do vetor A
    distribuicaoSuperficial = sum(amp)
    print("ρ =", distribuicaoSuperficial) # rho
    return distribuicaoSuperficial

# Distribuição superficial de carga para o caso inicializado (N == 20)
rho20 = distSuperf(A20)

ρ = 7.387195695632595e-07

# Resultado, gráfico da amplitude dos pulsos

def ampChart(amp):
    plt.title('Pulsos', fontsize=20)
    plt.xlabel('Elemento', fontsize=20) # Elemento = parte da placa dividida pelo segmento (N)
    plt.ylabel('Amplitude', fontsize=20) # Valor esperado
    plt.plot(amp, color='blue', linewidth='2')
    plt.show()

# Gráfico plotado para o caso inicializado (N == 20)
ampChart(A20)

# Gráfico 3D (distribuição de carga)

def ampChart3D(centros, amp, texto = 'Distribuição superficial de carga'):
    newR = np.array(centros)
    
    # Gráfico para z = 0
    # Primeira metade dos arrays
    X = newR[:len(newR)//2, 0]
    Y = newR[:len(newR)//2, 1]
    Z = np.array(amp[:len(amp)//2])

    figr = plt.figure(figsize=plt.figaspect(0.5));
    # ax = plt.axes(projection='3d');
    ax = figr.add_subplot(1,2,1, projection = '3d')
    c = X+Y+Z
    ax.scatter(X, Y, Z, c=c, cmap='viridis', edgecolor='none');
    ax.set_title('Distribuição de carga na placa com z = 0');
    # plt.show();

    # plt.figure();
    # ax = plt.axes(projection='3d');
    ax = figr.add_subplot(1,2,2, projection = '3d')
    ax.plot_trisurf(X, Y, Z, cmap='viridis', edgecolor='none');
    ax.set_title('Distribuição de carga na placa com z = 0');
    plt.show();

    print('z = d')
    print()

    # Gráfico para z = d
    # A outra metade dos arrays

    Xd = newR[len(newR)//2:len(newR), 0]
    Yd = newR[len(newR)//2:len(newR), 1]
    Zd = np.array(amp[len(amp)//2:len(amp)])

    # plt.figure();
    # ax = plt.axes(projection='3d');
    figr = plt.figure(figsize=plt.figaspect(0.5));
    ax = figr.add_subplot(1,2,1, projection = '3d')
    c= Xd+Yd+Zd
    ax.scatter(Xd, Yd, Zd, c=c, cmap='viridis', edgecolor='none');
    ax.set_title('Distribuição de carga na placa com z = d');
    # plt.show();

    # plt.figure();
    # ax = plt.axes(projection='3d');
    ax = figr.add_subplot(1,2,2, projection = '3d')
    ax.plot_trisurf(Xd, Yd, Zd, cmap='viridis', edgecolor='none');
    ax.set_title('Distribuição de carga na placa com z = d');
    plt.show();
    
# Gráfico 3D plotado para o caso inicializado (N == 20)
ampChart3D(matrizPosicao20, A20) 

z = d

Terceira questão

Com o aumento do número de segmentos (N), a carga tende a se concentrar na borda ("no canto") da placa

# Teste para vários N
listaN = [4, 8, 12, 20, 30, 40, 50, 60]
for k in listaN:
    delta = float(L)/k # lado de cada elemento
    qtdElementosK = int(3*((k**2)/2))
    Rk = coordenadasCentro(k)
    Zk = matrizImpedancia(qtdElementosK, Rk)
    Vk = matrizTensao(qtdElementosK)
    Ak = amplitude(Zk, Vk)
    print('Para N =', k)
    ampChart3D(Rk, Ak, k)

Para N = 4

z = d

Para N = 8

z = d

Para N = 12

z = d

Para N = 20

z = d

Para N = 30

z = d

Para N = 40

z = d

Para N = 50

z = d

Para N = 60

z = d

Quarta questão

Capacitância pode ser aproximada da seguinte forma:

$ C = \frac{Q}{V} $, onde $ V = V_{0} $ e $ Q = \Delta^{2} \cdot \sum_{i = N^{2}}^{2N^{2}} a_{i} $

O valor analítico pode ser obtido pela expressão abaixo:

$ C = \epsilon_{0} \frac{Área}{Distância} = \epsilon_{0} \frac{L^{2}}{d} $

# Cálculo da capacitância
def capacitancia(amp, n):
    # Não tem capacitância em z = 0 (tensão nula)
    return sum(amp[(n//2):n]) * delta**2/V0

capacitancias = []
alcance = range(5, 80, 5)
# 5, 10, ..., 70, 75
for k in alcance:
    delta = L/k # lado de cada elemento
    qtdElementosK = int(3*((k**2)/2))
    Rk = coordenadasCentro(k)
    Zk = matrizImpedancia(qtdElementosK, Rk)
    Vk = matrizTensao(qtdElementosK)
    Ak = amplitude(Zk, Vk)
    cap = capacitancia(Ak, qtdElementosK)
    capacitancias.append(cap)
    print('Capacitância para N =', k, 'é:', cap);

Capacitância para N = 5 é: -8.330586333689967e-13
Capacitância para N = 10 é: -2.989551097455612e-11
Capacitância para N = 15 é: 1.280470154347472e-10
Capacitância para N = 20 é: 5.464664783592663e-11
Capacitância para N = 25 é: 4.4924343508110533e-11
Capacitância para N = 30 é: 4.213417690321582e-11
Capacitância para N = 35 é: 4.0846181158309515e-11
Capacitância para N = 40 é: 4.0027663315589416e-11
Capacitância para N = 45 é: 3.983845572859667e-11
Capacitância para N = 50 é: 3.9292233466373106e-11
Capacitância para N = 55 é: 3.942084370155707e-11
Capacitância para N = 60 é: 3.891402877684328e-11
Capacitância para N = 65 é: 3.918459688830463e-11
Capacitância para N = 70 é: 3.867327666328806e-11
Capacitância para N = 75 é: 3.9026873227979175e-11

O objetivo é o erro (diferença entre as capacitâncias do analítico e computacional) ser próximo de 0.

plt.plot(np.array(alcance), np.array(capacitancias) - e0 * (3*(L**2)/4)/d)
plt.plot(np.array(alcance), np.zeros(len(alcance)))
plt.xlabel('N')
plt.ylabel('Capacitância')
plt.show()

# Valor analítico aproximado

print('Capacitância esperada =', e0 * (3*(L**2)/4)/d)

Capacitância esperada = 3.318750000000001e-11

Analisando o gráfico, é possível observar a convergência da capacitância com o valor analítico aproximado (esperado)